#include "nr.h"

void calcula_jacobiano(int nelements,
		       int nelectrodes,
		       int nnodes,
		       int CURR_PATTERN,
		       double **volt_pattern_matrix,
		       double corr,
		       int no_terra,
		       double **k_global,
		       double ***k_locais,
		       int *vec_electrodes,
		       double *distancias_centroides,
		       int nnp_local,
		       double *v_rho_0,
		       double ***H,
		       double *volt_calculado_diff,
		       int *VARIABLE_TABLE,
		       int *erro_saida)
{
  int i,m,n;
  int q_mais,q_menos;
  double *current;
  double *aux_line;
  double **k_global_inversa;
  double **Hj_diff;
  double **Hj_se;
  double *volt_calculado_SE,*volt_calculado_temp,aux;
  FILE *fp;
  
  k_global_inversa=dmatrix(1,nnodes,1,nnodes);
  aux_line=dvector(1,nnodes);
  
  for (m=1;m<=nnodes;m++){
    aux_line[m] = k_global[m][no_terra];
    k_global[no_terra][m]=0.0;
    k_global[m][no_terra]=0.0;
  }
  k_global[no_terra][no_terra]=1.0;

  invert_intel(nnodes,
	       k_global,
	       k_global_inversa,
	       erro_saida,
	       1);  /* 1: matriz simetrica   2:matriz nao simetrica*/

  /**********      restaura a matriz global          ********************/
  for (m=1;m<=nnodes;m++)
    k_global[no_terra][m] = k_global[m][no_terra] = aux_line[m];

  if((*erro_saida)!=0){
    free_dmatrix(k_global_inversa,1,nnodes,1,nnodes);
    free_dvector(aux_line,1,nnodes);
    (*erro_saida)=501;
    return;
  }

  volt_calculado_SE=dvector(1,nelectrodes);
  volt_calculado_temp=dvector(1,nelectrodes);
  Hj_diff=dmatrix(1,nelectrodes,1,nelements);
  (*H)=dmatrix(1,nelectrodes*nelectrodes,1,nelements);              VARIABLE_TABLE[20]=1;
    
  for (i=1;i<=nelectrodes;i++){
    Hj_se=dmatrix(1,nelectrodes,1,nelements);
    current=dvector(1,nnodes);
    /*printf("\t\t\t\tpadrao de corrente = %d \n",i);*/
    
    q_mais=i;
    q_menos=(q_mais + CURR_PATTERN)%nelectrodes+1;
    /*      q_terra=(q_mais + 1 + CURR_PATTERN)%nelectrodes+1;*/
    current[vec_electrodes[q_mais]]= corr;  
    current[vec_electrodes[q_menos]]= -corr;

    monta_sensibilidade (Hj_se,
			 k_global_inversa,
			 k_locais,
			 vec_electrodes,
			 current,
			 volt_calculado_SE, /* voltagens do problema direto SE*/
			 nnodes,
			 nelements,
			 nelectrodes,
			 nnp_local,
			 distancias_centroides, /**/
			 v_rho_0);

    /*******************   Calcula Hj_diff=volt_pattern_matrix*Hj_se ********************/

    SEtoDIF(volt_pattern_matrix,
	    Hj_se,
	    nelectrodes,
	    nelements,
	    Hj_diff);

    /*******************   calcula potencial calculado differencial ********************/

    multiplyV(nelectrodes,
	      nelectrodes,
	      &(volt_pattern_matrix[1][1]),  /*  Entrar com &(Matrix1[1][1]) ou &(Matrix1[0][0]) */
	      &(volt_calculado_SE[1]),  /*  Entrar com &(Vector1[1]) ou &(Vector1[0]) */
              NULL,
	      &(volt_calculado_temp[1]),  /*  Entrar com &(Vector2[1]) ou &(Vector2[0]) */
	      22);

    for (m=1;m<=nelectrodes;m++){
      for (n=1;n<=nelements;n++)
	(*H)[((nelectrodes)*(i-1))+m][n]=Hj_diff[m][n];
      volt_calculado_diff[nelectrodes*(i-1)+m]=volt_calculado_temp[m];
    }

    free_dmatrix(Hj_se,1,nelectrodes,1,nelements);
    free_dvector(current,1,nnodes);
  }
  
  print_dmatrix(nelectrodes*nelectrodes,nelements,(*H),"H_normal.txt");
  
  /*
  fp=fopen("idx_H.txt","r");
  fscanf(fp, "%lf", &aux);
  m=(int)aux;
  for (i=1;i<=m;i++){
    fscanf(fp, "%lf", &aux);
    n=(int)aux;
    for (m=1;m<=nelectrodes*nelectrodes;m++){
      (*H)[m][n]=0.0;
    }
  }
  fclose(fp);    
  */
  print_dmatrix(nelectrodes*nelectrodes,nelements,(*H),"H_com_zeros.txt");
  
      
  free_dvector(volt_calculado_SE,1,nelectrodes);
  free_dvector(volt_calculado_temp,1,nelectrodes);
  free_dmatrix(k_global_inversa,1,nnodes,1,nnodes);
  free_dvector(aux_line,1,nnodes);
  free_dmatrix(Hj_diff,1,nelectrodes,1,nelements);
}
